clc;clear;clearvars;
addpath('CEC2008\');
global initial_flag

% num_vari维度,max_evalutions最大评估次数,num_initial初始化点的个数,evalutions已经评估个数
num_vari = 100;
%max_evalutions = 5000 * num_vari;
max_evalutions = 500 * num_vari;
num_initial = 150;

% 搜索区间:f1[-100,100],f2[-100, 100],f3[-100,100],f4[-5,5],f5[-600,600],f6[-32,32],f7[-1,1],每一维搜索区间一样
upper_bound = [100,100,100,5,600,32,1];
lower_bound = [-100,-100,-100,-5,-600,-32,-1];

for func_num = 1:7

    % 使initial_flag==0满足此条件,换一个函数initial_flag重置为0
    initial_flag = 0;
    evalutions = 0;
    varphi = 0.1; %100维0,500维0.1和0.05,1000维0.15和0.1

    % 生成num_initial个种群,并获得其评估值
    sample_x = lhsdesign(num_initial, num_vari).*(upper_bound(func_num) - lower_bound(func_num)) + lower_bound(func_num).*ones(num_initial, num_vari);
    sample_y = benchmark_func(sample_x,func_num);
    evalutions = evalutions + size(sample_y,1);
    v = zeros(num_initial,num_vari);  %v初始化为0
    fmin = min(sample_y);
    arr_fmin = [];
    arr_fmin = [arr_fmin;fmin];

    while evalutions <= max_evalutions

        pb = sample_x;
        pby = sample_y;
        pbv = v;
        mean_x = mean(sample_x);

        next_samplex = [];
        next_sampley = [];
        next_v = [];

        % 更新pb
        while ~isempty(pb)
            k1 = randi(size(pb, 1));    %随机抽取三个点
            k2 = randi(size(pb, 1));
            while k2 == k1
                k2 = randi(size(pb, 1));
            end
            k3 = randi(size(pb, 1));
            while k3 == k2 || k3 == k1
                k3 = randi(size(pb, 1));
            end

            if pby(k1,1) < pby(k2,1)
                if pby(k3,1) < pby(k1,1)
                    add = k3;
                    update1 = k2;
                    update2 = k1;

                else
                    add = k1;
                    update1 = k2;
                    update2 = k3;
                end
            else
                if pby(k3,1) < pby(k2,1)
                    add = k3;
                    update1 = k2;
                    update2 = k1;
                else
                    add = k2;
                    update1 = k3;
                    update2 = k1;
                end
            end

            next_samplex = [next_samplex;pb(add,:)]; %add直接加入
            next_sampley = [next_sampley;pby(add,:)];
            next_v = [next_v;pbv(add,:)];
            pbv(update1,:) = rand(1, num_vari) .* pbv(update1,:) + rand(1, num_vari) .* (pb(add,:) - pb(update1,:)) + varphi .* rand(1, num_vari) .* (mean_x - pb(update1,:));    %update1需要更新
            xl = pb(update1,:) + pbv(update1,:);
            xl(xl > upper_bound(func_num)) = upper_bound(func_num);
            xl(xl < lower_bound(func_num)) = lower_bound(func_num);% 范围检查
            if evalutions <= max_evalutions
                yl = benchmark_func(xl,func_num);
                evalutions = evalutions + 1;
                if yl < fmin
                    fmin = yl;
                end
                arr_fmin = [arr_fmin;fmin];
                fprintf("func_num:%d evalutions:%d fmin: %.4f\n", func_num, evalutions, fmin);
                next_samplex = [next_samplex;xl];   %加到下一次迭代中
                next_sampley = [next_sampley;yl];
                next_v = [next_v;pbv(update1,:)];
            else
                break;
            end

            pbv(update2,:) = rand(1, num_vari) .* pbv(update2,:) + rand(1, num_vari) .* (pb(add,:) - pb(update2,:)) + varphi .* rand(1, num_vari) .* (mean_x - pb(update2,:));    %update2需要更新
            x2 = pb(update2,:) + pbv(update2,:);
            x2(x2 > upper_bound(func_num)) = upper_bound(func_num);
            x2(x2 < lower_bound(func_num)) = lower_bound(func_num);% 范围检查
            if evalutions <= max_evalutions
                y2 = benchmark_func(x2,func_num);
                evalutions = evalutions + 1;
                if y2 < fmin
                    fmin = y2;
                end
                arr_fmin = [arr_fmin;fmin];
                fprintf("func_num:%d evalutions:%d fmin: %.4f\n", func_num, evalutions, fmin);
                next_samplex = [next_samplex;x2];   %加到下一次迭代中
                next_sampley = [next_sampley;y2];
                next_v = [next_v;pbv(update2,:)];
            else
                break;
            end


            % 从pb中去掉k1,k2和k3,从后面删除,防止后面删除的数索引发生变化
            delete_i = [k1,k2,k3];
            delete_i = sort(delete_i,'descend');

            pb(delete_i(1), :) = [];
            pby(delete_i(1), :) = [];
            pbv(delete_i(1), :) = [];
            pb(delete_i(2), :) = [];
            pby(delete_i(2), :) = [];
            pbv(delete_i(2), :) = [];
            pb(delete_i(3), :) = [];
            pby(delete_i(3), :) = [];
            pbv(delete_i(3), :) = [];
        end

        % 更新sample_x,v和sample_y
        sample_x = next_samplex;
        sample_y = next_sampley;
        v = next_v;
    end
    plot(arr_fmin);
end